!the main program  (monte carlo simulation and parameter calculation)
subroutine exchange
use para
use mod_filename
use mpi
implicit none
!include 'mpif.h'


 character(len=8)::data1
 character(len=10)::time,zone
integer::num1,i,j,k
integer,dimension(1:8)::values
integer,dimension(1:8)::time1
integer,dimension(1:8)::time2
real::nmin,nsec,nhour
integer::isam

integer istatus(MPI_STATUS_SIZE)
integer,external::weight
double precision::timea,timeb,timec

!open(200,file='./time/time'//filename(my_rank)//'.dat')
!if(my_rank==0)open(200,file='./output/time.dat')
open(1000,file='./output/beta.dat')
open(1800,file='./output/timecount.dat')
open(500,file='./outdata/'//filename(my_rank)//'.bin',form='unformatted')
!open(1200,file='./output/betatun'//filename(my_rank)//'.dat')
open(3320,file='initialq.dat')

if(sam_flag==1)then
    if(my_rank==0)then
    do i=1,nbeta_sam
    open(5555+i,file='spin_sam'//filename(i)//'.bin',form='unformatted')
    open(6666+i,file='spin_sam'//filename(i)//'.dat')
    end do
    end if
end if

 !call MPI_INIT(ierr)
 !call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
 !call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
if(my_rank==0)then
 open(200,file='./output/time.dat')
 call date_and_time(data1,time,zone,values)
 write(200,*)data1,time,zone
 write(200,*)values
 write(200,'(10x)') 
 time1=values
end if 

 call readincar
 call latticedefine        !define the lattice number
 call neighbor             !define the neighbor lattice

 call config_init
 if(uflag==1)call initial_lamda
 call initial_energy
 call initial_mag
 call initial_count

 send_init(0)=tote/nlatt
 send_init(1)=totu/nlatt

! write(*, *) "TEST1" !20190305
 if(my_rank > 0) call MPI_SEND(send_init,2,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)
! write(*, *) "TEST2" !20190305
 
 if(my_rank==0)then
     write(3320,'(i5,3f16.5)')betanum(0),beta(betanum(0)),send_init(0),send_init(1)
 !write(*, *) "TEST3" !20190305
     do irank=1,nproc-1
 !write(*, *) "TEST4" !20190305
     call MPI_RECV(recv_init,2,MPI_DOUBLE_PRECISION,irank,itag,MPI_COMM_WORLD,istatus,ierr)
 !write(*, *) "TEST5" !20190305
     write(3320,'(i5,3f16.5)')betanum(irank),beta(betanum(irank)),recv_init(0),recv_init(1)
     end do
     
 end if
 
!equilibriation before exchange
 mybeta_num=betanum(my_rank)
 do iequi=1,nequi1
 call metropolis
if (my_rank==0 .and. mod(iequi, 500)==0) write(10000, "(A12,2X,I6)") "EQUILIBRIUM:", iequi 
 end do
!write(*, *) "test_1: ", tote  !20190305

!###################### added to test deltaE ####################

 !send_init(0)=tote/nlatt
 !send_init(1)=totu/nlatt

 !call MPI_SEND(send_init,2,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)

 !if(my_rank==0)then
 !    do irank=0,nproc-1
 !    call MPI_RECV(recv_init,2,MPI_DOUBLE_PRECISION,irank,itag,MPI_COMM_WORLD,istatus,ierr)
 !    write(3320,'(i5,3f16.5)')betanum(irank),beta(betanum(irank)),recv_init(0),recv_init(1)
 !    end do

 !end if

 !stop
!################################################################


!the recursion of the exchange process
  ischeme=0
 do irec=1,nrec
if(my_rank==0) write(10000,"(A10,2X,I3)") "RECURSION:",irec
!control the exchange times and exchange intervals    
    if(irec==1)then
    nexch=nexch1
    nswp=nswp1
    else
    nexch=nexch2
    nswp=nswp2
    end if

    exacpt(:)=1
    excount(:)=1
!equlibrium sweep


  do iexch=1,nexch
    !if(my_rank ==0 .and. mod(iexch, 100)==0) write(*, "(A9,2X,I6)") "EXCHANGE:", iexch
   if(iexch==1)timea=MPI_WTIME()  !count the time before equlibriation

    mybeta_num=betanum(my_rank)
    do iswp=1,nswp
    call metropolis
    if(my_rank ==0 .and. &
    & mod(iexch, 50)==0 .and. &
    & mod(iswp, 100)==0) write(10000, "(A9,2X,I6,A8,2X,I6)") &
    & "EXCHANGE:", iexch,"SWEEP:", iswp
    if(proflag==1) write(500) mybeta_num,tote,totm!, dot1e
    end do
!    write(500,'(2i8,f10.5)')mybeta_num,iexch,nacpt/ncount 

!################################################### output the configurations for MSD ###############################  

if(sam_flag==1)then 
if(mod(iexch,5)==1)then
!if(iexch>34500)then

do isam=1,nbeta_sam

if(mybeta_num==beta_sam(isam))then

if(my_rank > 0) call MPI_SEND(spin,nlatt*natom*3,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)

if(my_rank > 0) call MPI_SEND(udis,nlatt*nmod,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)

end if

if(my_rank==0)then
if (process(beta_sam(isam)) == 0) then
  write(5555+isam)iexch, 0 
else 
  call MPI_RECV(spin_out,nlatt*natom*3,MPI_DOUBLE_PRECISION,&
  &process(beta_sam(isam)),itag,MPI_COMM_WORLD,istatus,ierr)
  call MPI_RECV(udis_out,nlatt*nmod,MPI_DOUBLE_PRECISION,&
  &process(beta_sam(isam)),itag,MPI_COMM_WORLD,istatus,ierr)
  write(5555+isam)iexch,process(beta_sam(isam))
end if

!call MPI_RECV(spin_out,nlatt*natom*3,MPI_DOUBLE_PRECISION,process(beta_sam(isam)),itag,MPI_COMM_WORLD,istatus,ierr)
!call MPI_RECV(udis_out,nlatt*nmod,MPI_DOUBLE_PRECISION,process(beta_sam(isam)),itag,MPI_COMM_WORLD,istatus,ierr)
!  write(5555+isam)iexch,process(beta_sam(isam))
!   write(1100,*)'--------------'
   do i=1,nlatt
      do k=1,natom
      write(5555+isam)(spin_out(j,i,k),j=1,3)
      write(6666+isam,*)(spin_out(j,i,k),j=1,3)
      end do
   end do

   do i=1,nlatt
     write(5555+isam)(udis_out(i,j),j=1,nmod)
     write(6666+isam,*)(udis_out(i,j),j=1,nmod)
  end do

end if
   
end do

end if

end if

!####################################################################################################################

 if(exchange_flag==1)then !exchange_flag : whether to perform exchange

   if(iexch==1)timeb=MPI_WTIME()    !count time before exchange 

   if(ischeme>nscheme-1)ischeme=0
   ibeta=0
   do while(exbeta(ischeme,ibeta)>=0)
     tarbeta1=exbeta(ischeme,ibeta)
     tarbeta2=tarbeta1+1
     excount(tarbeta1)=excount(tarbeta1)+1

     !#### local block #####
     if(mybeta_num==tarbeta1)then

         if(my_rank/=process(tarbeta1))then
           write(*,*)'my_rank is',my_rank,'proess is',process(tarbeta1),'exchange is wrong'
           !call MPI_Finalize(ierr)
         end if

         nsend(0)=tote
        if(my_rank > 0) call MPI_SEND(nsend(0),1,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)
         !write(*,*)'first send is ok'
     end if
     
     if(my_rank==0)then
         !write(*,*)'the first recv is ok'
     end if
  
     if(mybeta_num==tarbeta2)then
         if(my_rank/=process(tarbeta2))then
         write(*,*)'my_rank is',my_rank,'proess is',process(tarbeta2),'exchange is wrong'
         !call MPI_Finalize(ierr)
         end if
         
         nsend(1)=tote
         !write(*, *) tote !20190305
        if(my_rank > 0) call MPI_SEND(nsend(0),1,MPI_DOUBLE_PRECISION,0,itag,MPI_COMM_WORLD,ierr)
         !write(*,*)'the second send is ok'
     end if

     if(my_rank==0)then
       if(process(tarbeta1) == 0) then
         nrecv(0) = nsend(0)
       else
         call MPI_RECV(nrecv(0),1,MPI_DOUBLE_PRECISION,process(tarbeta1),itag,MPI_COMM_WORLD,istatus,ierr)
       end if

       if(process(tarbeta2) == 0) then
         nrecv(1) = nsend(1)
       else
         call MPI_RECV(nrecv(1),1,MPI_DOUBLE_PRECISION,process(tarbeta2),itag,MPI_COMM_WORLD,istatus,ierr)
       end if

!     	 call MPI_RECV(nrecv(0),1,MPI_DOUBLE_PRECISION,process(tarbeta1),itag,MPI_COMM_WORLD,istatus,ierr)
!         call MPI_RECV(nrecv(1),1,MPI_DOUBLE_PRECISION,process(tarbeta2),itag,MPI_COMM_WORLD,istatus,ierr)
!         write(*,*)'the second recv is ok','nrecv=',nrecv(0),nrecv(1)
!         write(*,*)beta(tarbeta1),beta(tarbeta2)     
         exflag=weight(nrecv(1),nrecv(0),beta(tarbeta2),beta(tarbeta1))
!         write(*,*)'exflag is',exflag
         if(exflag==1)then
         betanum_temp=betanum(process(tarbeta1))
         betanum(process(tarbeta1))=betanum(process(tarbeta2))
         betanum(process(tarbeta2))=betanum_temp
         process_temp=process(tarbeta1)
         process(tarbeta1)=process(tarbeta2)
         process(tarbeta2)=process_temp
         exacpt(tarbeta1)=exacpt(tarbeta1)+1
         end if
     end if
     
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       call MPI_BCAST(betanum,nproc,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       call MPI_BCAST(process,nproc,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       call MPI_BCAST(exacpt,nproc,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) 
        
     !one pair of exchange is over
     ibeta=ibeta+1
        
     end do
!one scheme of exchange is over
     ischeme=ischeme+1

!     write(1200,*)iexch,betanum(my_rank)  ! output the beta tuneling rate 

     if(iexch==1)timec=MPI_WTIME()         !count time after one exchange
     
     if(iexch==1.and.my_rank==0)then
     write(1800,*)'the recursion is',irec
     write(1800,*)'the equlibration time is',timeb-timea,', the exchange time is',timec-timeb
     end if 

    end if   ! end if of exchange_flag

   end do

  if(exchange_flag==1)then
        
     if(my_rank==1)then
     write(*,*)'irec is ',irec,'the acpt distribution is as follows:'
     do i=0,nproc-1
     write(*,*)i,exacpt(i),excount(i)
     end do

     write(*,*)'process is as follows'
     do i=0,nproc-1
     write(*,*)i,process(i)
     end do
     end if

  end if
     
     if(proflag==0)then
        call recursion
        !equlibriation after recursion 

        mybeta_num=betanum(my_rank)
        do iequi=1,nequi2
        call metropolis
        end do
     end if
end do

!output the spin and udisp configration
        call spinout

 if(my_rank==0)then
      
 call date_and_time(data1,time,zone,values)
 write(200,*)data1,time,zone
 write(200,*)values
 time2=values


nmin=(time2(3)-time1(3))*24.0*60.0+(time2(5)-time1(5))*60.0+time2(6)-time1(6)+(time2(7)-time1(7))/60.0
nsec=nmin*60.0
nhour=nmin/60.0
write(200,'(10x)')
write(200,*)'total time used (min):',nmin
write(200,*)'total time used (hour):',nhour
write(200,*)'total time used (sec):',nsec


! close(100)
 close(200)

end if !end if my_rank==0

 close(500)
 close(1800)
 close(1000)

end subroutine

!the main program end 
!subroutine readincar
!use para
!implicit none
!open(2000,file='incar')
!read(2000,*)jexc
!read(2000,*)k31
!read(2000,*)k32
!read(2000,*)D3
!read(2000,*)D4
!if(hflag==1)then
!read(2000,*)hfield
!read(2000,*)magmom
!end if
!
!if(my_rank==0)then
!write(*,*)'k31,k32,D3,D4 are',k31,k32,D3,D4
!end if 
!!write(*,*)jexc,'jexc of',my_rank,' is ok'
!!write(*,*)j31,omig,'phonon of',my_rank,'is ok'
! close(2000)
!end
 

function weight(tote,tote2,beta,beta2) result(exflag)
implicit none
double precision,external::ran1
integer,save::weight_seed=-5000
double precision,intent(in)::tote,tote2,beta,beta2
integer::exflag
double precision::deltah
deltah=(beta-beta2)*(tote-tote2)
if(deltah>0.0.or.exp(deltah)>ran1(weight_seed))then
exflag=1
else
exflag=0
end if
end function
